Single hole dynamics in the t-J model on a square lattice 
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We present quantum Monte Carlo (QMC) simulations for a single hole in a t-J model from 
J = OAt to J — 4t on square lattices with up to 24 x 24 sites. The lower edge of the spectrum is 
directly extracted from the imaginary time Green's function. In agreement with earlier calculations, 
we find flat bands around (0, ±7r), (±7r, 0) and the minimum of the dispersion at (±7r/2, ±7r/2). For 
small J both self-consistent Born approximation and series expansions give a bandwidth for the lower 
edge of the spectrum in agreement with the simulations, whereas for J/t > 1, only series expansions 
agree quantitatively with our QMC results. This band corresponds to a coherent quasiparticle. This 
is shown by a finite size scaling of the quasiparticle weight Z(k) that leads to a finite result in the 
thermodynamic limit for the considered values of J/t. The spectral function A(k, to) is obtained 
from the imaginary time Green's function via the maximum entropy method. Resonances above the 
lowest edge of the spectrum are identified, whose J-dependence is quantitatively described by string 
excitations up to J/t — 2. 

PACS numbers: 71.10.Fd, 71.10.Pm, 75.10.Jm 



I. INTRODUCTION 

Since the pioneering work by Brinkman and RiceB 
(BR) the dynamics of a hole in an antiferromagnet re- 
mained as a recurring open problem in condensed mat- 
ter physics. After the discovery of high temperature 
superconductorsB and the suggestions by Andersonu on 
the possibility of a non-Fermi liquid state in those mate- 
rials, the question whether the quasiparticle weight of a 
hole vanishes due to the interaction with an antiferromag- 
netic background became central in the field of strongly 
correlated fermions. 

The BR treatment led to a fully incoherent spectrum in 
the so-called retraceable path approximation, for an anti- 
ferromagnetic Ising-like background, in the limit J™-^ 0. 
The retraceable path approximation is exact in oneu and 
in infinite dimensionsu but not in two dimensions since 
contributions of loops (Trugman pathaJ) may lead to 
a coherent propagation of the hole. Furthermore, for 
an Ising-like background it was shown within a Lanc- 
zos schemetj, that a finite quasiparticle weight is ob- 
tained. For the case of physical interest, namely with a 
Heisenberg spin-background, a large number of numer- 
ical methodsQ led to conflicting results. Whereas ex- 
act diagonalizations found large quasiparticle peaks at 
the lower edge of the spectrumO, quantum Monte Carlo 
(QMC) results were interpreted as leading to a vanish- 
ing quasiparticle weightEI. Since exact diagonalizations 
are possible only on very small lattices, finite size scaling 
cannot be performed reliably. On the other hand, QMC 
simulations suffered from the minus-sign problem, such 
that scaling was not possible with reasonable confidence. 
Further studies hasecLon the self-consistent Born approx- 
imation (SCBA)E9il3 gave a finite quasiparticle weight. 
However, since there fluctuations of the spin-background 
are only taken into account in the frame of a spin-wave 



approximation, the results obtained are not conclusive. 
Exact results for the-supersymmetric point J = 2t were 
obtained by SorellaE3, that give important benchmarks 
for any analytical or numerical method (see Sec. Ill B| ) , 
but unfortunately, they cannot be rigorously extended to 
the physical relevant parameter range J ~ OAt. 

Quite recently, the dynamics of a single hole in an an- 
tiferromagnetic background became experimentally ac- 
cessible by angle resolved photoemission spectroscopy 
(ARPES) in uadoped materials like Sr 2 Cu0 2 Cl 2 EiE3 and 
Ca 2 Cu0 2 Cl 2 lii The main features observed there are a 
minimum of the dispersion at k = (ir/2,ir/2) together 
with a vanishing of spectral weight beyond this point 
along the (1,1) direction. The obtained spectra show 
that the very flat portion around (ir, 0), that in optimally 
doped materials is almost degenerate with the bottom 
of the spectrum at (n/2, 7r/2)0, is shifted upwards (in 
a hole representation) by approximately 300 meV. This 
contradicts the single hole spectra found theoretically so 
far, where essentially the lower edge of the spectrum at 
k = (7r/2,7r/2) and (ir, 0) are almost degenerate, such 
that additional second andpt-hird nearest neighbor hop- 
ping terms were suggestedHL3, that lead to an agree- 
ment of the exact diagonalization results with experi- 
ments. Such terms were made recently responsible also 
for the vanishing of spectral weighLcbse to (ir/2, ir/2) by 
reducing the quasiparticle weightll3'E£l. 

In the following we present dynamical properties of 
a single hole in a two-dimensional t-J model on lat- 
tices with up to 24 x 24 sites in the parameter range 
0.4 < J/t < 4. Results were obtained with a new QMC 
algorithm, wherej-the spin background is simulated with 
a loop-algorithn£3 and the hole is exactly propagated 
for a given configuration of the spin background. The 
lower edge of the spectrum is obtained directly from the 
asymptotic form of the imaginary time Green's func- 



1 



tion. The resulting dispersion agrees with previous re- 
sults obtained within SCBA and series expansions 2 -] for 
J/t < 1, whereas for J ft > 1 only agreement with se- 
ries expansions is found. In particular, a flat disper- 
sion is obtained around k = (tt, 0) very close in value 
to the bottom of the baad|-at k — (tt/2,tt/2), in con- 
trast to the experimentsli-MJl The asymptotics of the 
imaginary time Green's function delivers also the quasi- 
particle weight for that band. Finite size scaling is pre- 
sented showing that Z(k) is finite for the parameter range 
considered, such that the lower edge of the spectrum 
corresponds to a coherent quasiparticle. Furthermore, 
our data are consistent with another exact predictionEj, 
namely that at the supersymmetric point and in the ther- 
modynamic limit, Z(Q)/Z{{)) = (2m) 2 , where Q = (tt, tt) 
is the antiferromagnetic wave vector and m is the stag- 
gered magnetization. The spectral function A(k,Lu) is 
calculated by analytic continuation with maximum en- 
tropy (MaxEnt)E 2 ]. Overall agreement is found with ex- 
act diagonalizations. At the supecsymmetric point, the 
delta function predicted by SorellaHl for the wave vector 
k = (0,0) is exactly reproduced. By extracting the con- 
tribution of the quasiparticle from the imaginary time 
Green's function, a resonance above the quasiparticle 
band is made evident, that together with the lower edge 
of the spectrum scales as J 2 / 3 , in agreement with the 
string pictures used to described the excitations for a 
hole in an antiferromagnetic Ising-background. Remark- 
ably, also the prefactors of the corresponding Airy func- 
tions are needed in order to properly describe the dis- 
tance between the resonance and the quasiparticle band. 

The paper is organized as follows. Section [n] describes 
the model, a canonical transformation that leads to a bi- 
linear form in spinless fermions interacting with S = h 
pseudospins, and the algorithm. Since the Hamiltonian 
for the transformed t-J model is bilinear in fermions (the 
holes) , their propagation can be calculated exactly given 
a pseudospin configuration. In Sec. [II the results are 
discussed. Section III A describes the lower edge of the 



spectrum and how it is obtained. In Sec. Ill B t he res ults 
for the quasiparticle weight are shown. Section III C de- 
scribes the spectral function A(k,uo) and the string exci- 



tations. Finally, the conclusions are given in Sec. IV 



II. THE MODEL AND THE ALGORITHM 

The t-J model is a suitable one to simulate the dynam- 
ics of a single hole in an antiferromagnet. On the one side, 
it can be obtained from the Hubbard model in the large 
coupling limit, which at half- filling leads to the Heiscn- 
berg antiferromagnet. On the other side, it is the rele- 
vant-one to simulate the cuprates, as shown by Zhang and 
RicecJ, and hence, to compare with experimentaLfr LLj. Its 
Hamiltonian is 
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neighbours only. In order to render Eq. (2.1) a bilin- 
ear form in fermionic operators, we perform a canonical 
transformationt3 
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where j h ± = (l±a i<z )/2 and er i)± = (a^ x ±iou y ) /2. The 
spinless fermion operators fulfill the canonical anticom- 
mutation relations {fj,fj} = Si,j, and (Ji^ a , a = x,y, or 
z, axe the Pauli matrices. The Hamiltonian becomes 
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<7j)/2, Ay = (1 - rii - 



and 



rii = flfi- The constraint to avoid doubly occupied 
states transforms to the conserved and holonomic con- 
straint Ei7*>— = 0- This constraint simply means, 
that a spinless fermion and a pseudospin J, are not al- 
lowed to sit on the same site. 

In order to obtain the dynamics of the hole, we calcu- 
late the one-particle Green's function for spin up, 

G(i-j,r) = -<Tc a (r)4 )T ) = -<T/t(r)/,.) (2.4) 



where T corresponds to the time ordering operator. In- 
serting complete sets of spin states the quantity above 
transforms as 
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TABLE I. Contributions for the propagation of the hole on one plaquette. The first column shows the weight for a 
propagation where the hole stays on the same site x, whereas in the second column the weight corresponds to the propagation 
to the adjacent site. The third column represents the spin background on the plaquette. 



Here mAr = f3, nAr — r, Art <C 1 and 
exp(— AtH(<ji, 02 )) is the evolution operator for the 
holes, given the spin configuration (01,02). I n the case 
of single hole dynamics, \v) is the vacuum state for holes, 
and P(a) is the probability distribution of a Heisenberg 
antiferromagnet for the configuration 0, where is a vec- 
tor containing all intermediate states (01, . . . a n , . . . , 01). 
The sum over spins is performed in a_very efficient way 
by using a world-line loop-algorithmEJ for a Heisenberg 
antiferromagnet with discretized imaginary time. In gen- 
eral we have At = 0.05, such that the extrapolation to 
At = leads to values of the observables within the sta- 
tistical error bars. The inverse temperature j3 is taken 
such that the energy is well converged (f3J > 20 for 
16 x 16 and (3 J = 30 for 24 x 24 sites), and there- 
fore, the data correspond to the ground-state. As the 
evolution operator for the holes is a bilinear form in 
the fermion operators, G(i,j,T,a) can be calculated ex- 
actly, in nontxast to a direct implementation in the loop 
algorithmcjE], where fermion paths are sampled stochas- 



tically. G(i, j, t, 0) contains a sum over all possible 
fermion paths between (i,0) and (j, r). The numerical 
effort to calculate G(i,j, t, 0) Vi,t scales as Nt, where 
N is the number of lattice points in space. Therefore, the 
present method is more efficient for large systems than 
e.g. projector algorithms for the Hubbard model, that 
scale with the system size cubed. 



With the representation of Eq. (2.3), the propagation 
of down spin electrons cannot be easily considered, since 
the operators 0^± cut world-lines. This is certainly not 
a problem for finite-size systems, where SU(2) symmetry 
is conserved. Since P(3) is the probability distribution 
for the quantum antiferromagnet, the algorithm does not 
suffer from sign problems on bipartite lattices and non- 
frustrating magnetic interactions in any dimension. 

We now address the explicit calculation of G(i,j, t, 0). 
In a first step, we introduce additional complete sets of 
single fermion states in the numerator of Eq.( |2.5[ ), such 
that G(i,j, t, 0) becomes 
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where the sum J"V runs 
ble intermediate one particle 
fermionic Hilbert space { | /; ) } . 

(fi p \e- A ^^ *-J\fi P ^)/(* P \e- 
only nonzero, when l p and l p ~\ belong to the same pla- 
quette. Therefore the entries of the matrices U (a p , P -i) 
are nonzero only at the positions which correspond to a 
plaquette in the checkerboard breakup. These entries are 
given in Table Q. As we are only interested in the Hilbert 
space of no double occupancy, we have to enforce the 



constraint at one single position of the propagation by 
projecting out the fermionic states which do not respect 
the constraint. We do so at r = corresponding to the 
first propagation. 

The one-particle spectral function 



A(k,u) = J2\(f,N-l\ % JO,N] 



(2.7) 
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is connected with the Green's function in imaginary time 
at T = 0, by the spectral theorem 



aa-.r) = I t fc; exp( Ttj) A(k,u). 



(2.8) 



spectral function A(k,ui) is discussed. Using the infor- 
mation from the lower edge of the spectrum, a resonance 
above the quasiparticle band is identified, that is very 
well described as a string excitation. 



Here | 0, TV) is the ground-state at half filling with energy 
Eq and | /, N — 1 ) are states in the N— 1 particle Hilbert 
space with energy E^ -1 . We perform the inversion of 

Eq. (pTsj), that due to the statistical errors of G(k,r) 
is an extremely ill-posed problem, by means of MaxEnt, 
where the A(k, ui) obtained is the one that maximizes the 
probability P(A\G), given the Green's function G(k,r). 
Correlations in the imaginary time data were taken into 
account by considering the covariance matrix. Details 
about MaxEnt can be found in the comprehensive review 
article by J.E. Gubernatis and M. JarrellE3. 

We would like to stress finally, that part of the dynam- 
ical data presented below were obtained without use of 
MaxEnt but directly extracted from the imaginary time 
Green's function. This is possible due to the high statis- 
tics and stability attainable with the present algorithm. 
The slowest decaying exponential, that corresponds to 
the excitation with lowest energy can be extracted sim- 
ply by fitting the tail of the Green's function at large 
values of r. This leads to the value of the excitation 



and its corresponding weight, as shown in Sees. HI A and 



[II B. Furthermore, in connection with MaxEnt, the next 



higher excitation can be obtained by subtracting the con- 
tribution from the lowest one from the Green's function. 
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FIG. 1. The energy of the lowest excitation is extracted 
from the imaginary time asymptotics of the Green's function, 
as indicated by the dotted line for a) J = 0.4t and b) J = 2t 
in a 16 x 16 lattice. 



This procedure is discussed in Sec. Ill C 



III. RESULTS 



A. The lower edge of the spectrum 



We concentrate in the following on three aspects of the 
dynamics of a single hole in a Heisen berg antiferromag- 
net. First we consider in Sec. Ill A the lower edge of 
the spectrum. This is a quantity that can be obtained 
by several other methods, including various Monte Carlo 
algorithms, such that the relative accuracy of each one 
and the region in parameter space, where each method 
gives best results, can be assessed. In our case, this quan- 
tity is obtained from the asymptotic behavior of the one- 
particle Green's function in imaginary time. However, 
not only the energy but also the weight of such an exci- 
tation can be extracted from the asymptotics, leadi ng to 
the quasiparticle weight, as discussed in Sec. Ill B. The 



present algorithm is up to now the only one capable of 
extracting this information for the t-J model free of ap- 
proximations on large lattices (in general up to 16 x 16 
and for J/t — 2 up to 24 x 24 sites). For small lattice 
sizes, the results can be compared with exact diagonal- 
izations, whereas for large systems only comparisons with 
approximate methods like SCBA can b e mad e. Finally, 
the whole spectrum is considered in Sec. HI C, where the 



The accuracy and stability of the data allow in our case 
to obtain the lower edge of the spectrum directly from the 
slope of the one-particle Green's function as a function 
of imaginary time r, for large values of r. Figure ^ shows 
the asymptotics in imaginary time for two values of the 
coupling constant, showing that the most accurate results 
are obtained, when J/t = 2. J/t = 0.4 is the smallest 
coupling, where such a procedure can be applied. In or- 
der to check the results obtained at the smallest coupling, 
we made additional calculations at Art = 0.2 (all other 
calculations are done at Art = 0.05), where larger values 
of rt can be reached. The resulting Green's functions 
are the same within the error bars, indicating a small At 
effect. 
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FIG. 2. Lower edge of the spectrum along the symmetry 
lines of the Brillouin zone for a) J/t — 0.4 and b) J/t — 2 in 
a 16 x 16 lattice. Comparisons are made with YMC (circles), 
GFMC for J/t = 0.4 (x), and series expansionO for J/t = 2 
(dotted line). 



suits are too high in energy, but still close to our numeri- 
cally exact ones. For values of J/t > 1, additional results 
from series expansions^ are available. At J/t = 2 we 
observe in general a very good agreement. Only around 
(tt, 0) we see that series expansions slightly underesti- 
mate the energy of the hole. The general features of 
the lower edge are not substantially modified when go- 
ing from J/t = 0.4 to J/t = 4. This is shown in Fig. 
||, where the only changes observed are an overall shift 
in energy with respect to the Heiscnbcrg antiferromagnet 
and a change in the bandwidth. The shift in energy can 
be followed by considering the dependence of e(n/2, tt/2) 
on J/t. This dependence is rather accurately described 
by e{J/t)/t = -3.28 + a^J/t) 2 / 3 , where a x is the first 
eigenvalue of th e dimensionless Airy equation (see Fig. 
|i"3| in Sec. [II C[ ). Such a scaling of the hole energies is 
found in th e tjJL m odel in the continuum limit for small 
values of J 2 cJii3ll3, when loops along the path of the hole 
are disregarded. In that case, the constant is —2\Jz — 1, 
where z is the coordination number. The resulting string 
picture gives an accurate description of the lowest exci- 
tations close to k = (tt/2, tt/2). As will be shown in Sec. 



Ill Q , also the next higher excitation can be described by 
the string picture. 



Figure |2j shows the lower edge of the spectrum for 
J/t = 0.4 and J/t = 2 in a 16 x 16 sites system. The 
energies are displayed with respect to the ground-state 
energy of the Heisenberg antiferromagnet. The results 
are compared with variational Monte Carlo (VMC)E3, 
Green's function Monte Carlo (GFMC)Ej, and series 
expansionsEU, whenever data is available. At J/t = 0.4 
(Fig. ||a), where our results are most affected by fluc- 
tuations, we observe good agreement with GFMC. The 
behavior of the statistical error is similar in both meth- 
ods, with larger fluctuations around k = (0, 0) and (tt, tt). 
Around k = (tt, 0) our results stow somewhat larger fluc- 
tuations. For J/t — 0.4 VMCeS, also appears to be very 
accurate concerning the lower edge. When its energies 
are compared to our calculations and the GFMC tech- 
nique, we find that their energies are within the error 
bars of the exact QMC calculations. At k = (0,0), 
the variational result is at the lower edge of the error 
bars of our calculation, and have the smallest statisti- 
cal error of all three approaches. At this specific fc-point 
both GFMC and our approach have large fluctuations 
before the state with lowest energy is clearly reached. As 
mentioned above, additional calculations with Art = 0.2 
were performed, in order to check the results obtained, 
without observing significant changes. 

Figure |^b shows that at J/t — 2, where our algorithm 
leads to much more accurate results, the variational re- 
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FIG. 3. Lower edge of the spectrum along the symmetry 
lines of the Brillouin zone for a) J/t = 0.6, b) J/t — 0.8, c) 
J/t = 1.2, and d) J/t = 4 in a 16 x 16 lattice. 
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FIG. 4. Bandwidth of the^ower edge as a function of J/t. 
e(0, 0) — e(7r/2,7r/2) in a 16 x 16 lattice compared with ex- 
act diagonalization (4x4 sites, open circles), GFMC (cross), 
VMC (open boxes), SCBA (dashed line), and series expan- 
sions (full line). 

Figure ^ shows the bandwidth obtained ip-pur simula- 
tions compared with exact diagonalizationstp, GFMCa, 
SCBABy, VMCE3 and series expansions^! For J/t < 
0.8 good agreement is found among all methods, whereas 
for larger values of J, only series expansions and VMC 
agree with our data. This, and the fact that the string 
picture gives a good representation of the lowest lying 
states, suggest that a perturbation expansion as per- 
formed in series expansions can be used to interpret the 
distinctive features of the lower edge. In particular the 



fiat band observed around k = (tt, 0) and the fact that the 
degeneracy between this poiat and k = (tt/2, 7r/2) sug- 
gested by some approaches!^ is clearly lifted, as shown 
by our simulation, are very well reproduced by series 
expansions. The flat bands can be well observed for 
all considered values of J/t, when considering the lower 
edge (Fig. and ||) and the complete spectral func- 
tion (Fig. [ll]). Our data clearly show for J/t > 0.6, that 
the neighboring points of k = (tt, 0) are generally slightly 
higher in energy. The band in this area does not seem 
to be completely flat, but it changes its curvature with 
local minima of the dispersion at the points (tt, 5) and 
(7T — 5, 0), when going along the (1,0) and (0,1) directions 
respectively, with the caveat that they are well defined 
beyond the error bars only for J/t > 1. In all the cases 
we find S rs 0.3-7T. This region with a very flat band spans 
an extremely large area in the Brillouin zone. 

A flat band on a similarly wide region in the Brillouin 
zone around k = (tt, 0) is also observed in photoemission 
spectroscopy of cuprates close to the Fermi-energy in the 
optimally doped compounds. As doping is reduced, that 
portion of the spectrum opens &-.pseudogap and weight 
is transferred to higher energieal3, until in the undoped 
materials, this portion is about 2 J (« 300meV) above 
the minimum at k = (7r/2, 7r/2)L3t3. The energy differ- 
ence between the points k — (tt/2, tt/2) and k = (tt, 0) 
is in our simulation about A = (0.25 ± 0.10)£ (sa J/2 
for J = OAt). The rather large error corresponds mainly 
to J/t < 1. No significant dependence on J/t can be 
observed in the whole range under consideration, in con- 
trast to the results from SCBA and series expansions. 
However, it could be that the J-dependence is masked in 
our case by large fluctuations, taking into account that 
the variations observed for this quantity by SCBA and se- 
ries expansions are much snialler than the one observed 
for the bandwidth. SCBAEil gives values ranging from 
0.17t (J/t = 1) to 0.12i (J/t = 4), that are smaller 
than the values we obtain. On the other hand, series 
expansionsEil obtain values between 0.15t at J/t = 1 and 
0.25£ at J/t = 2.5. The values obtained by series expan- 
sions are consistent with our results for large values of 
J/t. 



B. The quasiparticle weight 

The quasiparticle weight is the weight of the exponen- 
tial with the slowest decay, that is the exponential that 
determines the lower edge of the spectrum. This weight 
is 



Z(k)= lim G(-r,fc)exp[(eg-e )r] 



(3.1) 



In the following we focus on the thermodynamic limit 
of Z(k) for the wave vectors k = (tt, 0) and k = 
(tt/2, tt/2). Figures U and U show the finite-size scaling on 
these two points for J/t = 2 and J/t = 0.6 respectively. 
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For both k- and J-values, an appreciable quasiparticle 
weight is obtained, demonstrating that the lower edge of 
the spectrum describes the band of a coherent quasiparti- 
cle. The determination of the quasiparticle weight is only 
accurate for J/t > 0.6. Below that value, the quality of 
the data is less satisfactory and, for J/t = 0.4 the value 
presented can be taken only as an upper bound. The size 
dependence of Z(tt/2,tt/2) and Z(tt,0) is not very large 
and scales linearly with the inverse linear size-pf the sys- 
tem for J/t > 0.6, in agreement with SCBAHl. The size 
dependence at (tt/2, 7r/2) is systematically larger than at 
(tt, 0). The sizes considered are Lx L, with L = 16, 12, 8, 
and 4. At J/t = 2 we use additionally r-^pHi x 24 lat- 
tice. Values from exact diagonalizatiorHt 
eluded when available. 
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FIG. 5. a) Extrapolation 

of G(k, -r) = G(k, -r)exp[(e fc - e )r] for N = 8 x 8 and 
24 x 24 at J/t = 2. Finite-size scaling for b) k = (n/2, tt/2) 
and c) k = (tt, 0). 
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FIG. 7. Quasiparticle weight as a function of J/t for 
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exact diagonalizationEa in a 4 x 4 lattice is given by the tri- 
angle. 
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FIG. 8. Quasiparticle weight as a function of J/t for 
k = (7r/2,7r/2) (circles) and fc = (tt, 0) (x) in a 16 x 16 lat- 
tice. We compare our result with SCBA, where the dashed 
line corresponds to the quasiparticle weight for k = (n, 0) and 
the full line corresponds to k = (n/2,n/2). The data points 
were taken from Ref. 11. 
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from exact diagonalization resulta23'Eil. 



Figure [7| shows that the extrapolated quasiparticle 
weight increases with J/t both for k = (tt, 0) and k — 
(n/2, n/2). At J/t = 4 the quasiparticle reaches about 
80% of its maximal value. The changes of the quasipar- 
ticle weight with J/t are small when J/t > 1 and the 
slope becomes steeper for smaller values. Estimates of 
the quasiparticle weight were given both by VMCc3 and 
SCBAEH, the difference being rather small. The general 
trend is that VMC overestimates it at small J whereas 
SCBA overestimates it at large J. For definiteness we 
compare our results with SCBA for a 16 x 16 system 
in Fig. [|. We find a rather good agreement between 
both methods. As in our case Z(tt, 0) > Z(tt/2,tt/2) 
for all considered values of J/t. At small values of J 



7 



(0.01 < J/t < 0.5) SCBA finds a scaling of Z(tt/2, tt/2) = 
0.31 J 2/3 and Z(n,0) = 0.35 J 7 . For J/t > 1, the re- 
sults from SCBA overestimates the quasiparticle weight 
at the two considered fc-points, with an increasing devi- 
ation for larger values of J/t. Based on the quantitative 
agreement of SCBA with our results for small J, we can 
confidently conclude that the quasiparticle at k = (0, tt) 
and (7r/2,7r/2) should be finite for all values of J in the 
physically relevant region (i.e. J/i^O.l). 
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FIG. 9. Finite size scaling for the quasiparticle weight at 
Q — (n, 7r) for J/t — 2. The cross in the thermodynamic limit 
is (2m) 2 , m being the staggered magnetization. 



Z(k = 0) = 1/2, a requirement that is fulfilled by our 
simulation, where the Green's function at that particu- 
lar fc-point consists of a single exponential. In contrast 
to this, the estimate of SCBA is approximately 0.45 and 
that of VMC 0.32. Furthermore, Sorella showed that 
Z(Q)/Z(G) < (2m) 2 , where m 2 = S(Q)/N, S(Q) being 
the magnetic structure factor at the antifcrromagnetic 
wave vector. The equality is reached in the thermody- 
namic limit. Figure ^ shows the evolution with system 
size of Z(Q) together with results from exact diagonal- 
ization for a 4 x 4 system and (2m) 2 ~ 0.37 for L — > oo. 
Although large error bars show that the determination of 
Z(k) is less satisfactory for k = Q than at k = (tt/2, 7r/2), 
the data are coexistent with the exact result. It was fur- 
ther suggested^ 3 ! that if Z(k + Q)/Z(k) = (2m) 2 is satis- 
fied for k 7^ 0, a jump in the quasiparticle weight should 
be observed on crossing the border of the magnetic zone. 
Figure [l^ shows Z(k) along the symmetry directions in 
the Brillouin zone for a 24 x 24 system and J/t = 2. 
Our data do not show any sizable jump. Unfortunately, 
it is not possible to con sider arbitrarily long imaginary 
times since as Eq. (|3.l|) shows, the errors are amplified 
exponentially. Therefore, our results cannot be consid- 
ered as a proof of continuity. However, in view of the 
good agreement with the above mentioned exact results, 
we consider them as a convincing evidence. 



C. Spectral function and string excitations 




FIG. 10. Z(k) along the symmetry lines in the Brillouin 
zone for J/t = 2 in a 24 x 24 lattice. 

As mentioned in the introduction, there are exact 
results for the quasiparticle weight at the supersym- 
metric point in two dimcnsional3. On the one hand, 



The results discussed i n Sec . Ill A for the lower edge of 
the spectrum and in Sec. Ill B for the quasiparticle weight 
can be recognized in the spectral function (Fig. [ll]) ob- 
tained by using MaxEnt. For clarity, the maximum of 
each curve is normalized to 1 in the plots. The small 
numbers on the right hand side of the figures corre- 
spond to the maximal value of A{k 1 u) when the integral 

oo _ 

J dujA{k,uj) is properly normalized to tt/2. The lower 

— oo 

edge of the spectrum remains like in the previous section, 
but the accuracy of its location in A(k, to) is reduced by 
MaxEnt. The peaks around (0, it) and (7r/2, 7r/2) are gen- 
erally very sharp, in agreement with the fact that a finite 
quasiparticle weight was found in Sec. [II B. A trans- 



fer of weight from high to low energies can be observed, 
when J/t is increased, consistent wi th th e increase in the 
quasiparticle weight (Fig. in Sec. | IIIB| ). 

When compared to the ID caseEj, it is seen that the 
high energy excitations in the 2D case are extremely 
broad. The total bandwidth remains essentially constant 
as a function of J in contrast to the ID case, where it 
scales as 4i + J. 
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FIG. 11. Spectral function for a 16 x 16 system and 
a) J/t = 0.4, b) 0.6, c) 0.8, d) 1.2, e) 2 and f) 4. The verti- 
cal lines indicate resonances abov e the quasi par ticle peak at 
k = (tt/2, tt/2) as obtained in Sec. [ill C| , Fig.]!^ 



For values of the coupling in the range J/t < 2 we ob- 
serve satellite peaks in the region around k = (tt/2, tt/2) 
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(Fig. [fl]) next to the lowest energy peak which is ex- 
tremely sharp and corresponds to a quasiparticle. The 
S-peak can not be handled satisfactorily by MaxEnt. As 
can be seen by comparison of Figs. || and 11, MaxEnt 
gives some weight at energies lower than the band edge. 
This additional weight has to be balanced in some way, 
such that this error propagates to the other side of the 
(5-peak. Small peaks in the vicinity of the (S-peak can 
therefore not be resolved. In order to resolve structures 
close to the quasiparticle peak, we subtract the exponen- 
tial corresponding to the lowest energy (see Fig. |l^). The 
thus modified Green's function can now be used as input 
of MaxEnt . 
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FIG. 12. Original Green's function and Green's function 
with subtraction (lower curve) for a 16 x 16 system and 
J/t = 2 at k = (7r/2,7r/2). 



(I.I) 



(|,0) 



FIG. 13. Original spectrum (full line) and after subtrac- 
tion (dashed line) along the line (n/2,k y ). In the direction 
toward (ir/2,0) the lowest resonance approaches the position 
of the quasiparticle peak and merges with it, whereas toward 
(7r/2,7r) the distance stays approximately constant. In the 
second case, the main effect is a broadening of the resonances. 
Shown is a 16 x 16 system with J/t = 0.6. 
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FIG. 14. Green's function and resulting spectral function 
before and after subtraction of the lowest exponential, where 
the lowest exponential is in the error bars of section [II A 



The results in (b) correspond (from bottom to top) the orig- 
inal spectral function, the spectral function used in Fig. ha, 
the result when subtracting the lowest possible exponential, 
and to the highest possible exponential. The error bar on the 
exponential is 0.03i, the error on its weight is 0.025. 

Before proceeding to the results, let us remark that, 
when the MaxEnt results obtained with the modified 
Green's function, i.e. after the subtraction of the lowest 
exponential, are viewed closely, on occasions, an addi- 
tional peak appears at the bottom of the spectrum (this 
effect can be seen e.g. for (7r/2,7r/2) in Fig. [l3]). To 
exclude, that this peak corresponds to a real physical 
effect, we take several modified Green's functions, that 
are consistent with the exponential of the lowest peak, 
within the statistical error. Therefore, we take the low- 
est and the highest expone ntial, t hat are consistent with 
the results obtained in Sec. Ill A, and use them as input 
of MaxEnt. As can be seen in Fig. 0, the new peak 
that appears below the low energy peak of the original 
function, and hence is artificial, is only observed in two 
cases with varying position, whereas the two other peaks 
can be always observed, no matter which exponential is 
subtracted (always within the statistical errors). The po- 
sition of these high energy peaks is not changed by the 
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different subtractions, only the width is affected. In all 
cases discussed, a small shift of these structures can be 
observed with respect to the ones in the spectrum with- 
out the subtraction. However, the positions assumed by 
these structures after the subtraction is not affected by 
the different subtractions within the values allowed by 
the statistical errors. We conclude that the initial small 
shift is due to the inability of MaxEnt to concentrate the 
weight of the delta- function of the quasiparticle peak to 
a single energy value. 



the position of the maxima of the existing peaks at small 
J/t compared to the full spectral function (except the 
low energy peak, that disappeared). We can further ob- 
serve, that the satellite peak next to the low energy peak 
can now be seen for all values of J/t < 2. One should 
notice, that no additional weight has been produced at 
high energies, but the normalization has changed (again 
the maximal value is normalized to 1, not the area of the 
spectral function). 
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FIG. 15. Change of the spectral function when subtracting 
the first excitation for a) (tt/2, tt/2) and b) (tt, 0). The energy 
is shifted by an arbitrary amount £, in order to display the 
spectra for different values of J/t in the same energy range. 
The dotted line is the original result, the full line gives the 
modified one. 

The result of the procedure described above is shown 
in Fig. [l5| For (ir/2, tt/2) there are only little changes of 




FIG. 16. The first three excitations at k = (tt/2, tt/2). At 
J/t = 0.4 only two peaks were resolved. The lines represent 
the solutions obtained by solving the linear string potential 
for the hole in the t — J z model. 

At (tt/2, 7t/2) (Fig. |l^(a)) the resolution of the second- 
lowest excitation is quite clear, when applying the above 
method, whereas at (tt,0) (Fig. |l5|(b)) the results are 
either not accurate enough, or the corresponding exci- 
tation is weaker. For J/t — 1.2 the resolution is not 
good enough to separate the two resonances at (w, 0) . 
Generally the excitations at higher energies at (tt, 0) are 
broader than at (tt/2, tt/2), so that the positions of the 
maxima are not as well defined. .Similar structures, were 
observed in exact diagonalizationuLI and in SCBAE3, and 
were ascribed to string excitations. 

When the string picture is valid, as it is expected in the 
t-J z model the hole is confined by a linear potential, [lead-, 
ing to (k-independent) eigenvalues of the energyE3t3'Lil 
given by 



Ejt = -2\/3 + a„(J 2 /i) 2/3 , 



(3.2) 



where a—, are the eigenvalues of a dimcnsionless Airy 
equationllj. The first three eigenvalues are given by 
a n = 2.33,4.08,5.52. In Fig. f^the results for the first 
three excitations are given for k = (tt/2, tt/2), and are 
compared to the predictions from SCBA. The error bars 
on the second and third peak are obtained as the width 
of the MaxEnt peak at half inte nsity, the error bars of 
the first peak are taken as in Sec. [II A| . We find, that for 
J/t < 2 the lowest peak can be accurately described by 
e (7r/2, tt/2) = -E H - 3.28t + 2.33( J/tflH, where E H is 
the Heisenberg energy per site, and the second peak by 
ei(7r/2,7r/2) = -E H -3.28t + 4:.08(J/t) 2 /H. The value of 
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3.28£ + Eh is the result obtained from SCBA0, whereas 
the prefactors of (J/t) 2 / 3 are exactly the values of the 
dimensionless Airy function, implying that the first two 
peaks behave (within our error bars) exactly as it is ex- 
pected by the string picture. In contrast to this, a fit 
from SCBA for the first three excitations in the t-J model 
for values of J/t < 0.4 results in a n = 2.16,5.46,7.81, 
also with the exponent 2 / 3li3, leading to a clear disagree- 
ment with our data. The third peak that can be resolved 
cannot be explained by the string picture, since its dis- 
tance to the lower band edge is independent of J and 
has a value of about At. The existence of a string ex- 
citation is not restricted to (tt/2,tt/2), but it can also 
be observed between (tt/2 7 it/2) and (n/2, 3ir/A). This is 
demonstrated for the value J/t = 0.6 (see Fig. |l3f) . 

The results above lead to the conclusion that the lowest 
excitations can be well described by the string picture. 
However, it should be kept in mind that the string picture 
originates in the Ising limit for J/t <C 1, and that it is 
based on the continuum limit, that seems far away from 
our case with strings of lengths between two and a maxi- 
mum of five lattice points, that correspond to the first two 
string excitations. Moreover, the string picture predicts 
a band without dispersion, that is clearly not the case in 
our simulations. A way to reconcile this paradoxical sit- 
uation is given by the very good quantitative agreement 
between QMC and series expansionsEil for the dispersion 
of the quasiparticle and its bandwidth for a fairly large 
range in J. As shown by the expansion around the Ising 
limit, a coherent motion of the hole is made possible af- 
ter the creation of strings due to hopping processes, by 
appropriate spin-flips, the shortest string being of length 
two. The lowest order contribution appears in third or- 
der, where the points (ir/2,ir/2) and (tt, 0) are degener- 
ate. Fourth and higher order processes remove this de- 
generacy, giving rise to a band that agrees qualitatively 
very well with the one obtained in QMC. Therefore on 
top of the coherent motion determined by J± , string like 
excitations are possible and related to J z and t. Such a 
possibility was already proposed by Beran, Poilblanc and 
Laughlin (BPL)Eil on the basis of exact diagonalizations 
on small systems and is confirmed unambiguously by our 
simulations on large systems. 

At J/t ^ 2 the excitations fall below the values pre- 
dicted by the string picture. In those regions the string 
picture is no longer valid, as the relaxation of the dis- 
turbed spin bonds is faster than the motion of the hole. 



IV. CONCLUSIONS 



A new QMC algorithm was presented that allows a 
rather accurate determination of the single hole dynam- 
ics in a two dimensional Heisenberg S — 1/2 antiferro- 
magnet. The main advantages of this algorithm are the 
combination of the loop algorithm for the update of the 
spins and the exact evolution of the hole for a given spin 



configuration. Due to the diverging correlation length at 
zero temperature, large autocorrelation times should be 
expected for algorithms with local updates, a problem 
that is avoided here by the global update of the spins. 
On the other hand, the exact evolution of the hole for 
a given spin background avoids further statistical errors 
that would be introduced if the hole is updated stochasti- 
cally, as in recently proposed approachesEljtZl. In fact, the 
accuracy achieved allows for a determination of several 
dynamical quantities on large lattice sizes, leading to the 
possibility of a finite size scaling of e.g. the quasiparticle 
weight. 



First (Sec. Ill A) we discussed the lower edge of the 
spectrum that is obtained directly from the asymptotics 
in imaginary time of the Green's function. This quantity 
is accessible to different techniques, that are however, 
with the exception of GFMC, either restricted to small 
lattice sizes or approximate. The comparison shows that 
very accurate results are given by series expansions over 
a large range of parameters, supporting thus the inter- 
pretation of the relevant physical processes for the co- 
herent motion of the hole in the frame of a pertubative 
expansion around the Ising limit. This picture is further 
enfor ced by our study of the quasip articl e weight (Sec. 
[II B) and the spectral function (Sec. [II C| ). In Sec. Ill B 



it was shown, that indeed the lower edge of the spectrum 
describes the coherent propagation of a hole with finite 
quasiparticle weight. This is the case for all the param- 
eter range studied, and due to the good agreement with 
SCBA especially for small values of J, one can conclude 
that this coherent propagation takes place for essentially 
J > 0. Furthermore, by considering structu res ne xt to 
the lowest peak in the spectral function (Sec. [II C), it is 
seen that the lowest excitations around the wave vector 
k = (n/2, 7r/2) are very well described by the levels of 
strings usually discussed for the t-J z model, giving fur- 
ther support to the pertubative picture, where the hole 
creates strings during its motion through the lattice, that 
are healed by exchange processes, leading thus to coher- 
ence. In fact, the strings for the first two levels, that 
agree quantitatively with our simulations correspond to 
a length of two and five lattice sites. Strings of length two 
are the dominant contributions in series expansions for 
the dispersion of the quasiparticle. Moreover, our find- 
ings showing the existence of string resonances above the 
quasiaarticle pole lend support to a picture developed by 
BPD 3 -!!, where the composite nature of the quasiparticle 
is advanced. In previous exact diagonalization studies, 
the existence of such resonances, that were first observed 
in 4 x A latticesL 2 ] were not clea rly id entified on larger 
latticesEl. We have shown in Sec. [II C that they can be 
quantitatively identified with string excitations. How- 
ever, they are visible only in a rather narrow region along 
the line k x ~ 7r/2, ir/2 < k y < 37r/4, such that in small 
lattices with up to 26 sites, these features can be very 
much affected by boundary effects. Following BPL, the 
quasiparticle can be viewed as a light holon attached to 
a spinon by a confining potential, the one that gives rise 
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to the spectrum of string excitations. 

A comparison with experiments0~Ei fails due to the 
small gap between the lowest peak at k = (ir/2,ir/2) 
and the flat band around k = (tt, 0). ft was suggested 
by several authorsE-3ll3llj that this shift might be ob- 
tained introducing hopping terms to second and third 
nearest neighbop. Furthermore, it was found in exact 
diagonalizationsEj that such extra terms lead to a no- 
ticeable reduction of the quasiparticle weight. Since ex- 
act diagonalizations with second and third nearest neigh- 
bors in lattices with 18 and 26 sites suffer considerably 
under finite size effects, a discussion of the influence of 
longer range hopping on the quasiparticle weight must 
be carried out in much larger lattices. Such studies are 
presently under way. 
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